#---------------------------------------------------------------------------------------------------------# 
#---------------------------------------------------------------------------------------------------------# 
# ------------------------------------        0. Max Joosten         ------------------------------------ #
#-------------------------------------   Last updated: 08.09.2023    -------------------------------------# 
#---------------------------------------------------------------------------------------------------------# 

# Appendix Tables A5; B1 - B5; C2; F1 - F3; H2 - H4

# 0. Load relevant packages
library(haven)       # For reading foreign data (read_stata)
library(tidyverse)   # Collection of useful packages (mostly dplyr)
library(kableExtra)  # For exporting tables
library(stargazer)   # For printing regression tables
library(sandwich)    # For adjusting standard errors

# 0. Set working directory
setwd("C:/Users/MAJOOST/OneDrive - UvA/phd/PhD project/002. MJ WS paper/01_data")
df <- read_stata("df_complete.dta")

# ------------------------------------          Table A5         ------------------------------------ #
# Voter preferences and Party positions per party, year and area
# Part I: Party 
tab_party <- df %>% group_by(party) %>%
  summarise(meani1 = mean(ipref_1, na.rm = T),
            meani2 = mean(ipref_2, na.rm = T),
            meani3 = mean(ipref_3, na.rm = T),
            meanpp = mean(pp, na.rm = T)*10)

# Two decimals
tab_party <- round(tab_party,2)

# Rename columns
colnames(tab_party) <- c("Party", "Lower", "Middle", "Higher", "Party position")

# Rename rows
tab_party <- tab_party %>% mutate(Party = case_when(Party == 1 ~ "CDA",
                                                    Party == 2 ~ "PvdA",
                                                    Party == 3 ~ "SP",
                                                    Party == 4 ~ "VVD",
                                                    Party == 5 ~ "GL",
                                                    Party == 6 ~ "CU",
                                                    Party == 7 ~ "D66",
                                                    Party == 8 ~ "SGP",
                                                    Party == 9 ~ "PVV"))

# Export table
tab_party %>%
  kbl(caption = "Table A5_1: Voter preferences by party") %>%
  kable_classic(full_width = F, html_font = "Times New Roman") %>% 
  save_kable("../05_figures/tablea5_1.html")

# Part II: Year 
tab_year <- df %>% group_by(year) %>%
  summarise(meani1 = mean(ipref_1, na.rm = T),
            meani2 = mean(ipref_2, na.rm = T),
            meani3 = mean(ipref_3, na.rm = T),
            meanpp = mean(pp, na.rm = T)*10)

# Add column names (for saliinequ variable) - Look at codebook for labels!
colnames(tab_year) <- c("Year", "Lower", "Middle", "Higher", "Party position")

tab_year <- round(tab_year, 2)

tab_year %>%
  kbl(caption = "Table A5_2: Voter preferences by year") %>%
  kable_classic(full_width = F, html_font = "Times New Roman") %>% 
  save_kable("../05_figures/tablea5_2.html")

# Part III: Area
tab_area <- df %>% group_by(area) %>%
  summarise(meani1 = mean(ipref_1, na.rm = T),
            meani2 = mean(ipref_2, na.rm = T),
            meani3 = mean(ipref_3, na.rm = T),
            meanpp = mean(pp, na.rm = T)*10)

# Add column names (for saliinequ variable) - Look at codebook for labels!
colnames(tab_area) <- c("Area", "Lower", "Middle", "Higher", "Party position")
tab_area <- round(tab_area, 2)
tab_area <- tab_area %>% mutate(Area = case_when(   Area == 2 ~ "Environment",
                                                    Area == 3 ~ "Development",
                                                    Area == 4 ~ "Education",
                                                    Area == 5 ~ "Transport",
                                                    Area == 6 ~ "Health care",
                                                    Area == 7 ~ "Crime",
                                                    Area == 8 ~ "Defence",
                                                    Area == 9 ~ "Social security"))


tab_area %>%
  kbl(caption = "Table A5_3: Voter preferences by area") %>%
  kable_classic(full_width = F, html_font = "Times New Roman") %>% 
  save_kable("../05_figures/tablea5_3.html")

## Appendix B: Different model specifications
# ------------------------------------          Table B1         ------------------------------------ #
# Adaptation in two surveys following the elections
mb1 <- lm(lead_ipref_1 ~ pp + ipref_1 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_1, data = df)
mb2 <- lm(lead_ipref_2 ~ pp + ipref_2 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_2, data = df)
mb3 <- lm(lead_ipref_3 ~ pp + ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_3, data = df)

# Adjust standard errors
robust_se_mb1 <- sqrt(diag(vcovHC(mb1, type = "HC1")))
robust_se_mb2 <- sqrt(diag(vcovHC(mb2, type = "HC1")))
robust_se_mb3 <- sqrt(diag(vcovHC(mb3, type = "HC1")))

# Output
stargazer(mb1, mb2, mb3,  
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),  
          out = "../05_figures/appendix_table_B1.html",
          dep.var.labels = c("Lower income", "Middle income", "Higher income"),
          se = list(robust_se_mb1, robust_se_mb2, robust_se_mb3),
          header=FALSE,
          title = "Table B1: Adaptation at t + 2",
          no.space = TRUE,
          column.sep.width = "3pt",
          font.size = "small",
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Party positions", "Lagged preference (at t + 2)"),
          add.lines=list(c("Fixed effects", "Yes", "Yes", "Yes")))

# ------------------------------------          Table B2         ------------------------------------ #
# Influence, different specifications
mb4 <- lm(pp ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + as.factor(area) + spendchange, weights = iwt123, data = df)
mb5 <- lm(pp ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + as.factor(year) + spendchange, weights = iwt123, data = df)
mb6 <- lm(pp ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + as.factor(party) + spendchange, weights = iwt123, data = df)
mb7 <-lm(pp ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, data = df)
mb8 <- lm(pp ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + l_pp + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt123, data = df)

# Adjust standard errors
robust_se_mb4 <- sqrt(diag(vcovHC(mb4, type = "HC1")))
robust_se_mb5 <- sqrt(diag(vcovHC(mb5, type = "HC1")))
robust_se_mb6 <- sqrt(diag(vcovHC(mb6, type = "HC1")))
robust_se_mb7 <- sqrt(diag(vcovHC(mb7, type = "HC1")))
robust_se_mb8 <- sqrt(diag(vcovHC(mb8, type = "HC1")))

# Output
stargazer(mb4, 
          mb5,
          mb6,
          mb7,
          mb8,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),
          out = "../05_figures/appendix_table_B2.html",
          dep.var.labels = c("Party positions"),
          se = list(robust_se_mb4,
                    robust_se_mb5,
                    robust_se_mb6,
                    robust_se_mb7,
                    robust_se_mb8),
          header=FALSE, 
          title = "Table B2: Linear regression models of parties' spending positions with different model specifications",
          no.space = TRUE, 
          column.sep.width = "3pt", 
          font.size = "small", 
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Lower income voter preferences", "Middle income voter preferences", "Higher income voter preferences", "Lagged party position", "Change in spending since previous election"),
          add.lines=list(c("Specification", "Area FE", "Year FE", "Party FE", "No weights", "Lagged PP")))

# ------------------------------------     Table B3 - B5         ----------------------------------- #
# Adaptation, different specifications
mb9  <- lm(ipref_1 ~ pp + l_ipref_1 + spendchange, weights = iwt_1, data = df)
mb10 <- lm(ipref_2 ~ pp + l_ipref_2 + spendchange, weights = iwt_2, data = df)
mb11 <- lm(ipref_3 ~ pp + l_ipref_3 + spendchange, weights = iwt_3, data = df)

mb12 <- lm(ipref_1 ~ pp + l_ipref_1 + as.factor(area) + spendchange, weights = iwt_1, data = df)
mb13 <- lm(ipref_2 ~ pp + l_ipref_2 + as.factor(area) + spendchange, weights = iwt_2, data = df)
mb14 <- lm(ipref_3 ~ pp + l_ipref_3 + as.factor(area) + spendchange, weights = iwt_3, data = df)

mb15 <- lm(ipref_1 ~ pp + l_ipref_1 + as.factor(year) + spendchange, weights = iwt_1, data = df)
mb16 <- lm(ipref_2 ~ pp + l_ipref_2 + as.factor(year) + spendchange, weights = iwt_2, data = df)
mb17 <- lm(ipref_3 ~ pp + l_ipref_3 + as.factor(year) + spendchange, weights = iwt_3, data = df)

mb18 <- lm(ipref_1 ~ pp + l_ipref_1 + as.factor(party) + spendchange, weights = iwt_1, data = df)
mb19 <- lm(ipref_2 ~ pp + l_ipref_2 + as.factor(party) + spendchange, weights = iwt_2, data = df)
mb20 <- lm(ipref_3 ~ pp + l_ipref_3 + as.factor(party) + spendchange, weights = iwt_3, data = df)

mb21 <- lm(ipref_1 ~ pp + l_ipref_1 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, data = df)
mb22 <- lm(ipref_2 ~ pp + l_ipref_2 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, data = df)
mb23 <- lm(ipref_3 ~ pp + l_ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, data = df)

# Adjust standard errors
robust_se_mb9  <- sqrt(diag(vcovHC(mb9, type = "HC1")))
robust_se_mb10 <- sqrt(diag(vcovHC(mb10, type = "HC1")))
robust_se_mb11 <- sqrt(diag(vcovHC(mb11, type = "HC1")))

robust_se_mb12 <- sqrt(diag(vcovHC(mb12, type = "HC1")))
robust_se_mb13 <- sqrt(diag(vcovHC(mb13, type = "HC1")))
robust_se_mb14 <- sqrt(diag(vcovHC(mb14, type = "HC1")))

robust_se_mb15 <- sqrt(diag(vcovHC(mb15, type = "HC1")))
robust_se_mb16 <- sqrt(diag(vcovHC(mb16, type = "HC1")))
robust_se_mb17 <- sqrt(diag(vcovHC(mb17, type = "HC1")))

robust_se_mb18 <- sqrt(diag(vcovHC(mb18, type = "HC1")))
robust_se_mb19 <- sqrt(diag(vcovHC(mb19, type = "HC1")))
robust_se_mb20 <- sqrt(diag(vcovHC(mb20, type = "HC1")))

robust_se_mb21 <- sqrt(diag(vcovHC(mb21, type = "HC1")))
robust_se_mb22 <- sqrt(diag(vcovHC(mb22, type = "HC1")))
robust_se_mb23 <- sqrt(diag(vcovHC(mb23, type = "HC1")))

# B3
stargazer(mb9,
          mb12,
          mb15,
          mb18,
          mb21,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),  
          out = "../05_figures/appendix_table_B3.html",
          dep.var.labels = c("Lower income", "Middle income", "Higher income"),
          se = list(robust_se_mb9,
                    robust_se_mb12, 
                    robust_se_mb15, 
                    robust_se_mb18, 
                    robust_se_mb21),
          header=FALSE,
          title = "Table B3: Linear regression models of the lower income spending preferences with different model specifications",
          no.space = TRUE,
          column.sep.width = "3pt",
          covariate.labels=c("Party positions", "Lagged lower income", "Lagged middle income", "Lagged high income", "Change in spending since previous election"),
          add.lines=list(c("Specification", "No FE", "Area FE", "Year FE", "Party FE", "No weights")))

# B4
stargazer(mb10,
          mb13,
          mb16,
          mb19,
          mb22,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),  
          out = "../05_figures/appendix_table_B4.html",
          dep.var.labels = c("Lower income", "Middle income", "Higher income"),
          se = list(robust_se_mb10,
                    robust_se_mb13, 
                    robust_se_mb16, 
                    robust_se_mb19, 
                    robust_se_mb22),
          header=FALSE,
          title = "Table B4: Linear regression models of the middle income spending preferences with different model specifications",
          no.space = TRUE,
          column.sep.width = "3pt",
          covariate.labels=c("Party positions", "Lagged lower income", "Lagged middle income", "Lagged high income", "Change in spending since previous election"),
          add.lines=list(c("Specification", "No FE", "Area FE", "Year FE", "Party FE", "No weights")))

# B5
stargazer(mb11,
          mb14,
          mb17,
          mb20,
          mb23,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),  
          out = "../05_figures/appendix_table_B5.html",
          dep.var.labels = c("Lower income", "Middle income", "Higher income"),
          se = list(robust_se_mb11,
                    robust_se_mb14, 
                    robust_se_mb17, 
                    robust_se_mb20, 
                    robust_se_mb23),
          header=FALSE,
          title = "Table B5: Linear regression models of the higher income spending preferences with different model specifications",
          no.space = TRUE,
          column.sep.width = "3pt",
          covariate.labels=c("Party positions", "Lagged lower income", "Lagged middle income", "Lagged high income", "Change in spending since previous election"),
          add.lines=list(c("Specification", "No FE", "Area FE", "Year FE", "Party FE", "No weights")))

## Appendix C2: Education analysis
mc1 <- lm(pp ~ l_epref_1 + l_epref_2 + l_epref_3 + as.factor(area) + as.factor(party) + as.factor(year) + spendchange, weights = ewt123, data = df)
mc2 <- lm(pp ~ epref_1 + epref_2 + epref_3 + as.factor(area) + as.factor(party) + as.factor(year) + spendchange, weights = ewt123, data = df)
mc3 <- lm(epref_1 ~ pp + l_epref_1 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = ewt_1, data = df)
mc4 <- lm(epref_2 ~ pp + l_epref_2 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = ewt_2, data = df)
mc5 <- lm(epref_3 ~ pp + l_epref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = ewt_3, data = df)

# Adjust standard errors
robust_se_mc1 <- sqrt(diag(vcovHC(mc1, type = "HC1")))
robust_se_mc2 <- sqrt(diag(vcovHC(mc2, type = "HC1")))
robust_se_mc3 <- sqrt(diag(vcovHC(mc3, type = "HC1")))
robust_se_mc4 <- sqrt(diag(vcovHC(mc4, type = "HC1")))
robust_se_mc5 <- sqrt(diag(vcovHC(mc5, type = "HC1")))

# Output
stargazer(mc1,
          mc2,
          mc3,
          mc4,
          mc5,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),
          out = "../05_figures/appendix_table_C2.html",
          dep.var.labels = c("Party positions", "Lower educated", "Middle educated", "Higher educated"),
          se = list(robust_se_mc1,
                    robust_se_mc2,
                    robust_se_mc3,
                    robust_se_mc4,
                    robust_se_mc5),
          header=FALSE, 
          title = "Table C2: Education analyses",
          no.space = TRUE, 
          column.sep.width = "3pt", 
          font.size = "small", 
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Party positions", 
                             "Lower educated voter preferences (t - 1)", "Middle educated voter preferences (t - 1)", "Higher educated voter preferences (t - 1)", 
                             "Lower educated voter preferences (t + 1)", "Middle educated voter preferences (t + 1)", "Higher educated voter preferences (t + 1)", 
                             "Change in spending since previous election"),
          add.lines=list(c("Fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes")))

# Cleaning global environment
rm(list=setdiff(ls(), c("df")))

# Appendix F: Validating the CPB party position measure
# Correlations with MP / CHES
df_mpchestable <- read_stata("mpchestable_may23.dta")

# Table F1
cor_tab <- df_mpchestable %>% select(!c(party, year)) %>%
  cor(use = "complete.obs") %>%
  round(2) %>%
  as.data.frame()

setwd("C:/Users/MAJOOST/OneDrive - UvA/phd/PhD project/002. MJ WS paper/01_data/Appendix")
write.xlsx(cor_tab,"appendix_table_F1.xlsx", row.names = T)

# Analyses with MP / CHES
df_mpches <- read_stata("mpches_may23.dta")

# Influence
mf1 <- lm(pp_mp ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt123, data = df_mpches)
mf2 <- lm(pp_ches ~ l_ipref_1 + l_ipref_2 + l_ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt123, data = df_mpches)

# Adjust standard errors
robust_se_mf1 <- sqrt(diag(vcovHC(mf1, type = "HC1")))
robust_se_mf2 <- sqrt(diag(vcovHC(mf2, type = "HC1")))

# Output
stargazer(mf1,
          mf2,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),
          out = "../05_figures/appendix_table_F2.html",
          dep.var.labels = c("Party positions (CMP)", "Party positions (CHES)"),
          se = list(robust_se_mf1,
                    robust_se_mf2),
          header=FALSE, 
          title = "Table F2: Linear regression models of parties' spending positions using MPD and CHES party positions",
          no.space = TRUE, 
          column.sep.width = "3pt", 
          font.size = "small", 
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Lower income voter preferences", "Middle income voter preferences", "Higher income voter preferences", "Change in spending since previous election"),
          add.lines=list(c("Fixed effects", "Yes", "Yes", "Yes", "Yes")))

rm(list=setdiff(ls(), c("df_mpches")))

# Table F3
mf3 <- lm(ipref_1 ~ pp_mp + l_ipref_1 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_1, data = df_mpches)
mf4 <- lm(ipref_2 ~ pp_mp + l_ipref_2 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_2, data = df_mpches)
mf5 <- lm(ipref_3 ~ pp_mp + l_ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_3, data = df_mpches)
mf6 <- lm(ipref_1 ~ pp_ches + l_ipref_1 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_1, data = df_mpches)
mf7 <- lm(ipref_2 ~ pp_ches + l_ipref_2 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_2, data = df_mpches)
mf8 <- lm(ipref_3 ~ pp_ches + l_ipref_3 + as.factor(area) + as.factor(year) + as.factor(party) + spendchange, weights = iwt_3, data = df_mpches)

# Adjust standard errors
r_se_mf3 <- sqrt(diag(vcovHC(mf3, type = "HC1")))
r_se_mf4 <- sqrt(diag(vcovHC(mf4, type = "HC1")))
r_se_mf5 <- sqrt(diag(vcovHC(mf5, type = "HC1")))
r_se_mf6 <- sqrt(diag(vcovHC(mf6, type = "HC1")))
r_se_mf7 <- sqrt(diag(vcovHC(mf7, type = "HC1")))
r_se_mf8 <- sqrt(diag(vcovHC(mf8, type = "HC1")))

# Output
stargazer(mf3,
          mf4, 
          mf5,
          mf6,
          mf7, 
          mf8,
          type = "text", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),  
          #out = "../05_figures/table_appendix_F3.html",
          dep.var.labels = c("Low (MPD)", "Middle (MPD)", "High (MPD)",
                             "Low (CHES)", "Middle (CHES)", "High (CHES)"),
          se = list(r_se_mf3, 
                    r_se_mf4, 
                    r_se_mf5,
                    r_se_mf6,
                    r_se_mf7,
                    r_se_mf8),
          header=FALSE,
          title = "Linear regression models of the public's spending preferences using MPD and CHES party positions",
          no.space = TRUE,
          column.sep.width = "3pt",
          #font.size = "tiny",
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Party positions (MPD)", "Party positions (MPD)", "Lagged lower income", "Lagged middle income", "Lagged high income", "Change in spending since previous election"),
          add.lines=list(c("Fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

rm(list=setdiff(ls(), c("df_mpches")))

# Appendix H: Household income
# Table H2: Household income (correlation)
df_hhinc <- read_stata("cob_essentials_hhincmay23.dta") %>% select(wave, party, area, l_hhipref_1, l_hhipref_2, l_hhipref_3)
df <- read_stata("cob_essentials_may23.dta") %>% filter(year %in% 2017) %>% select(wave, party, area, l_ipref_1, l_ipref_2, l_ipref_3)

df_hhcor <- df %>% left_join(df_hhinc,
                        by = c("wave", "party", "area")) %>%
  select(starts_with("l_"))

cor(df_hhcor)

# Table H3: Household income (regression, influence)
df_hhinc <- read_stata("cob_essentials_hhincmay23.dta")

# 1. Influence income (Stata code for exact same results: reg pp l_ipref_1 l_ipref_2 l_ipref_3 i.area i.year i.party spendchange [aweight=iwt123], vce(rob)
mh1 <- lm(pp ~ l_hhipref_1 + as.factor(area) + as.factor(party), weights = hhiwt_1, data = df_hhinc)
mh2 <- lm(pp ~ l_hhipref_2 + as.factor(area) + as.factor(party), weights = hhiwt_2, data = df_hhinc)
mh3 <- lm(pp ~ l_hhipref_3 + as.factor(area) + as.factor(party), weights = hhiwt_3, data = df_hhinc)
mh4 <- lm(pp ~ l_hhipref_1 + l_hhipref_2 + l_hhipref_3 + as.factor(area) + as.factor(party), weights = hhiwt123, data = df_hhinc)

# Adjust standard errors
robust_se_mh1 <- sqrt(diag(vcovHC(mh1, type = "HC1")))
robust_se_mh2 <- sqrt(diag(vcovHC(mh2, type = "HC1")))
robust_se_mh3 <- sqrt(diag(vcovHC(mh3, type = "HC1")))
robust_se_mh4 <- sqrt(diag(vcovHC(mh4, type = "HC1")))

# Output
stargazer(mh1, 
          mh2, 
          mh3, 
          mh4, 
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),
          out = "../05_figures/appendix_table_H3.html",
          dep.var.labels = c("Party positions"),
          se = list(
            robust_se_mh1, 
            robust_se_mh2, 
            robust_se_mh3, 
            robust_se_mh4),
          header=FALSE, 
          title = "Table H3: Linear regression models of parties' spending positions using household income",
          no.space = TRUE, 
          column.sep.width = "3pt", 
          font.size = "small", 
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Lower income voter preferences (HH income)", "Middle income voter preferences (HH income)", "Higher income voter preferences (HH income)"),
          add.lines=list(c("Fixed effects (Party / Area)", "Yes", "Yes", "Yes", "Yes")))

# Table H3: Household income (regression, adaptation)
mh5 <- lm(hhipref_1 ~ pp + l_hhipref_1 + as.factor(area) + as.factor(party), weights = hhiwt_1, data = df_hhinc)
mh6 <- lm(hhipref_2 ~ pp + l_hhipref_2 + as.factor(area) + as.factor(party), weights = hhiwt_2, data = df_hhinc)
mh7 <- lm(hhipref_3 ~ pp + l_hhipref_3 + as.factor(area) + as.factor(party), weights = hhiwt_3, data = df_hhinc)

# Adjust standard errors
robust_se_mh5 <- sqrt(diag(vcovHC(mh5, type = "HC1")))
robust_se_mh6 <- sqrt(diag(vcovHC(mh6, type = "HC1")))
robust_se_mh7 <- sqrt(diag(vcovHC(mh7, type = "HC1")))

# Output
stargazer(mh5, mh6, mh7,
          type = "html", 
          omit = c("party", "year", "area"),
          keep.stat=c("n", "adj.rsq"),  
          out = "../05_figures/appendix_table_H4.html",
          dep.var.labels = c("Lower income (HH income)", "Middle income (HH income)", "Higher income (HH income)"),
          se = list(robust_se_mh5, robust_se_mh6, robust_se_mh7),
          header=FALSE,
          title = "Table H4: Linear regression models of the public's spending preferences using household income",
          no.space = TRUE,
          column.sep.width = "3pt",
          #font.size = "tiny",
          notes = "All models include fixed effects for years, parties and policy areas.",
          covariate.labels=c("Party positions", "Lagged lower income (HH income)", "Lagged middle income (HH income)", "Lagged high income (HH income)"),
          add.lines=list(c("Fixed effects", "Yes", "Yes", "Yes")))
